# LOADING ####
# packages ####
library(ggmap)
library(igraph)
library(ggplot2)
library(ggrepel)
library(dplyr)
library(reshape2)
library(xtable)

# data ####
#Archaeological societies dataset:
archsoc.df <- read.csv2("archaeological-societies-cths.csv", stringsAsFactors = F)
archsoc.df <- archsoc.df[ ! is.na(archsoc.df$creation),]

#Anthropological institutions dataset:
inst1953.df <- read.csv2("anthropo-institutions1953.csv", stringsAsFactors = F)

# France base map ####
basemap <- map_data(map = "france")   
basemap <- ggplot(basemap) + theme_void()
basemap <- basemap + geom_polygon(aes(x = long, y = lat, group = group), 
                          colour="gray40", fill="white", size = .2 )
basemap <- basemap + coord_fixed(xlim = c(-4.5, 7.9), ylim = c(51,42), ratio = 1.5)  

# cities geocoding:
cities.df <- data.frame(location = unique(archsoc.df$location), stringsAsFactors = F)
cities.df <- cbind(cities.df, geocode(cities.df$location, source = "dsk"))

archsoc.df <- merge(archsoc.df, cities.df, by="location", all.x = T)
archsoc.df <- archsoc.df[order(archsoc.df$creation),]

archsoc.df$scope <- factor(archsoc.df$scope,
                          levels=c("multidisciplinary", "archaeology", "prehistory"))

# GROWTH OF ARCHAEOLOGICAL SOCIETES ####
# Foundations by year: ####
archsoc.by.year.df <- as.data.frame(table(archsoc.df$creation), stringsAsFactors = F)
colnames(archsoc.by.year.df) <- c("year", "n")
archsoc.by.year.df$year <- as.numeric(archsoc.by.year.df$year)
archsoc.by.year.df$year <- floor(archsoc.by.year.df$year / 10) * 10
archsoc.by.year.df <- as.data.frame(table(archsoc.by.year.df$year))
colnames(archsoc.by.year.df) <- c("year", "n")
print(xtable(archsoc.by.year.df, align="lll"),   include.rownames=F)

# Foundations by year and scope: ####
creations.by.scope <- archsoc.df %>% group_by(creation, scope) %>% summarise(n = n() )

creations.by.scope <- creations.by.scope %>% 
                    group_by(scope) %>% mutate(count = cumsum(n))

ggplot(creations.by.scope, aes(x = creation, y = count)) +
    theme_light(base_size=9, base_family = "Courier") +
    geom_point(stat="identity", colour="gray40", size=.8) +
    geom_path(size=.3) +
    scale_x_continuous("years", breaks = seq(1825, 2000, 25)) +
    ggtitle("Growth of archaeological societies in France",
            subtitle="Cumulative number by year and scope (n=104)") +
    theme(plot.title = element_text(hjust = 0.5)) +
    facet_wrap(~scope, nrow=3) +
    theme(plot.title = element_text(hjust = 0.5, margin=unit(c(.5,0, .2,0), "cm"),
                                    size= 10, face="bold"),
          plot.subtitle = element_text(hjust = 0.5, size= 8),  
          plot.margin=unit(c(0,.2, 0, 0), "cm")
          ) 
ggsave("archsoc-time-scope.pdf",  width=11, height = 9, units="cm", scale= 1)
ggsave("archsoc-time-scope.tiff", width=11, height = 9, units="cm", scale= 1, dpi = 300)
 
# GEOGRAPHY OF ARCHAEOLOGICAL SOCIETES ####

# Map of archaeological societies, 19th and 20th centuries ####
archsoc.df$dec <- floor(archsoc.df$creation / 10) * 10
archsoc.df$dec <- factor(archsoc.df$dec)
archsoc.df <- archsoc.df[order(archsoc.df$dec, decreasing = F ),]

basemap + 
    theme_void(base_size=7, base_family = "Courier") +
    geom_jitter(data = archsoc.df, aes(x=lon, y=lat,  fill=dec, shape=scope ), 
                alpha=.8, size=2, width=.07, height=.07 ) + 
    scale_fill_manual("Foundation decade:", breaks = seq(1820, 2000, 20),
            values = grey.colors(length(unique(archsoc.df$dec)),  start = 0, end = 1) ) +
    scale_shape_manual("Scope:", values=c(21,22,23) ) +
    ggtitle("Growth of archaeological societies\n in France (1824-2004)",
            subtitle="Decade of foundation and location (n=104)") +
    theme(legend.position="bottom",
          legend.box = "vertical",
          legend.box.just= "left",
          legend.key = element_rect(size =2, fill="white"),
          legend.key.size = unit(3.6, 'mm'),
          plot.title = element_text(hjust = 0.5, margin=unit(c(.5,0, .2,0), "cm"),
                                    size= 10, face="bold" ),
          plot.subtitle = element_text(hjust = 0.5, size= 8),
          plot.margin = unit(c(-.2, 0, 0, 0), "cm"),
          legend.title = element_text(size=7,  face="bold"),
          legend.direction = "horizontal",
          legend.margin = margin(t= -.2, r=0, b=0, l=0, unit="cm"),
          legend.box.margin = margin(t= -.8, r=0, b=0, l=0, unit="cm")
    ) +
    guides(fill = guide_legend(order = 2, nrow=1, title.position="top",
                               override.aes = list(shape = 21 )),
           shape = guide_legend(order = 1, title.position="top") ) 

ggsave("archsoc-map.pdf",  width=11, height = 14, units="cm", scale= 1)
ggsave("archsoc-map2.tiff",  width=11, height = 14, units="cm", scale= 1, dpi = 300)

# Map of archaeological societies in 1953 ####
archsoc1953.df <- archsoc.df[which(archsoc.df$creation <= 1953),]
archsoc1953.df$dec <- floor(archsoc1953.df$creation / 10) * 10
archsoc1953.df$dec <- factor(archsoc1953.df$dec)
archsoc1953.df <- archsoc1953.df[order(archsoc1953.df$dec, decreasing = F ),]

basemap + 
    theme_void(base_size=7, base_family = "Courier") +
    geom_jitter(data = archsoc1953.df, aes(x=lon, y=lat, shape = scope,  fill=dec), 
                alpha=.8, size= 2, 
                width=.06, height=.06 ) + 
    scale_fill_manual("Foundation decade:", breaks=seq(1820, 1950, 20),
                      values=grey.colors(13,  start = 0, end = 1) ) +
    scale_shape_manual("Scope:", values=c(21,22,23) ) +
    ggtitle("Archaeological societies in France in 1953",
            subtitle="Decade of foundation and scope (n=40)") +
    theme(legend.position="bottom",
          legend.box = "vertical",
          legend.box.just= "left",
          legend.key = element_rect(size =2, fill="white"),
          legend.key.size = unit(3.6, 'mm'),
          plot.title = element_text(hjust = 0.5, margin=unit(c(.5,0, .2,0), "cm"),
                                    size= 10, face="bold" ),
          plot.subtitle = element_text(hjust = 0.5, size= 8),
          plot.margin=unit(c(-1, 0, -.5, 0), "cm"),
          legend.title=element_text(size=7,  face="bold"),
          legend.margin=margin(t= 0, r=0.1, b=0, l=0, unit="cm"),
          legend.box.margin=margin(t= -1.2, r=0, b=0, l=0, unit="cm")
    ) +
    guides(fill = guide_legend(order = 2, nrow=1, title.position="top",
                           override.aes = list(shape = 21 )),
          shape = guide_legend(order = 1, title.position="top") ) 

ggsave("archsoc1953-map.pdf",  width=11, height = 14, units="cm", scale= 1)
ggsave("archsoc1953-map.tiff",  width=11, height = 14, units="cm", scale= 1, dpi = 300)

# ANTHROPOLOGICAL INSTITUTIONS IN 1953 ####

# Selection of the revelant institutions:
inst1953.df <- inst1953.df[inst1953.df$prehist.archaeo.related  == "yes",]

cities.df <- data.frame(location = unique(inst1953.df$location), stringsAsFactors = F)
cities.df <- cbind(cities.df, geocode(cities.df$location, source = "dsk"))

inst1953.df2 <- inst1953.df %>% group_by(location, type) %>% summarise(n = n() )
inst1953.df2  <- merge(inst1953.df2, cities.df, by="location", all.x = T)

# map of the prehistoric archaeology organizations in 1953 ####
basemap + 
    theme_void(base_size=7, base_family = "Courier") +
    geom_point(data = inst1953.df2, aes(x=lon, y=lat, size=n, shape = type),
                alpha=.6,  colour = "black", fill = "gray60") + 
    scale_shape_manual("Organization types:",  values = c(21,22,24))+ 
    # scale_size_area("Count:", max_size = 12, breaks = c(1,3, 6) ) +
    scale_size("Count:", range = c(4, 10), breaks = c(1,3, 6) ) +
    geom_label_repel(data = cities.df, aes(x=lon, y=lat, label=location), 
                     nudge_x = 0, nudge_y = 0.4, size=2,  label.padding = unit(0.1, "lines"),
                     min.segment.length=unit(.99, "lines"), label.size= 0 ) +
    ggtitle("Organizations related to prehistoric\narchaeology in France in 1953",
            subtitle="Typology and organization count by place (n=24)") +
    theme(legend.position="bottom",
          legend.box = "horizontal",
          legend.box.just= "top",
          legend.direction = "horizontal",
          legend.key = element_rect(size = 2, fill="white"),
          legend.key.size = unit(4, 'mm'),
          plot.title = element_text(hjust = 0.5, size= 10, face="bold" ),
          plot.subtitle = element_text(hjust = 0.5, size= 8),
          legend.margin = margin(t= -.2, r=0, b=0, l=0, unit="cm"),
          legend.title = element_text(size = 7,  face="bold"),
          plot.margin=unit(c(0.1 ,0, 0, 0), "cm"),
          legend.box.margin=margin(t= -.7, r=0, b=0, l=0, unit="cm")
          ) +
    guides(shape = guide_legend(order = 1, nrow=3, title.position="top",
                                override.aes = list(size = 3)),
           size = guide_legend(nrow=1, title.position="top"),
           type = guide_legend(keywidth = .5, keyheight = .5)) 

ggsave("institutions1953-map.pdf",  width=10, height = 13, units="cm", scale= 1)
ggsave("institutions1953-map.tiff",  width=10, height = 13, units="cm", scale= 1)

#  make bipartite graph (institution-actors): ####
inst.actors.df <- melt(inst1953.df[, c(2, 8:26)], id.vars="name2")
inst.actors.df <- inst.actors.df[,-2]
colnames(inst.actors.df) <- c("institution", "actor")
inst.actors.df <- inst.actors.df[ inst.actors.df$actor != "", ]
inst.actors.df$actor <- tolower(inst.actors.df$actor)
inst.actors.df$actor <- sub("(.*_.)([a-z]*)", "\\1", inst.actors.df$actor)
inst.actors.df$actor <-   sub("(.*_.*-.).*", "\\1", inst.actors.df$actor)

inst.actors.g <- graph.data.frame(inst.actors.df, directed=T)
V(inst.actors.g)$type <- "institution"
V(inst.actors.g)[ V(inst.actors.g)$name  %in% inst.actors.df$actor ]$type <- "actors"

# study of the relations between institutions ####
# — projection ####
inst.mat <- bibcoupling(inst.actors.g)
inst.g <- graph.adjacency(inst.mat, mode="undirected")
inst.g <- delete.vertices(inst.g, V(inst.g)$name %in%  inst.actors.df$actor )
E(inst.g)$weight <- 1
inst.simple.g <- simplify(inst.g,  edge.attr.comb="sum")

# — compute weighted vertex degree:####
strength.res <- sort(strength(inst.simple.g), decreasing = T)
xtable(as.data.frame(strength.res), digits= 0)


# — descriptive stats ####
bi.graphs.l <- list("Actors network"=actors.g, "Institutions network"=inst.g)
g.stats <- data.frame(
    Nodes = sapply(bi.graphs.l, gorder),
    Edges = sapply(bi.graphs.l, ecount),
    Components = sapply(bi.graphs.l, function(x) clusters(x)$no),
    "Degree centralization" = sapply(bi.graphs.l,
                                     function(x) round(centr_degree(x, mode="all")$centralization, 2)),
    "Betweeness centralization" = sapply(bi.graphs.l,
                                         function(x) round(centr_betw(x)$centralization, 2))
)

g.stats <- t(g.stats)
# xtable(g.stats)

# study of the relations between actors ####
# — projection  ####
actors.mat <- cocitation(inst.actors.g)
actors.g <- graph.adjacency(actors.mat, mode="undirected")
actors.g <- delete.vertices(actors.g, V(actors.g)$name  %in% inst.actors.df$institution)
E(actors.g)$weight <- 1
actors.simple.g <- simplify(actors.g, edge.attr.comb="sum")

# — compute betweeneness centrality ####
V(actors.simple.g)$size <- betweenness(actors.simple.g, directed= F, normalized=T) * 100

# — identification of the articulations points ####
articulation_points(actors.simple.g)


# — plot graph with ggplot####
gg <- delete_vertices(actors.g, degree(actors.g) == 0)
gg <-  decompose.graph(gg)[[1]]
V(gg)$size <- betweenness(gg, directed= F, normalized=T) * 100
    
v <- data.frame(name=V(gg)$name, betw=V(gg)$size, layout_nicely(gg))
v[order(v$betw),]

colnames(e) <- c("x", "y")
e <- data.frame(get.edgelist(gg))
colnames(e) <- c("name1", "name2")
e <- merge(e, v[,-2], by.x="name1", by.y="name", sort = F)
colnames(e) <- c("name1", "name2", "x1", "y1")
e <- merge(e, v[,-2], by.x="name2", by.y="name", sort = F)
colnames(e) <- c("name1", "name2", "x1", "y1", "x2", "y2")

v[v$betw == 0,]$ name <- NA
v[v$betw < sort(v$betw, decreasing = T)[11],]$ name <- NA

svg("actors.graph.svg")
ggplot() +
    theme_nothing() +
    geom_segment(data = e, aes(x= x1, y= y1, xend =x2, yend= y2), color="grey") +
    geom_point(data=v, aes(x = X1, y = X2), size=1) +
    geom_label(data= v, aes(x = X1, y = X2,  label = name), size =2) +
    ggtitle("")   
dev.off()

